#include "dyn_specialFunctions_besselj.h"

BesselJ::BesselJ(double _order)
{
    setOrder(_order);
}

double BesselJ::operator ()(const double x)
{
    double nu = fabs(order);

    if (order < 0.0)
    {
       return boost::math::cyl_bessel_j(nu,x)*cos(arma::datum::pi*nu)
               -boost::math::cyl_neumann(nu,x)*sin(arma::datum::pi*nu);
    }
    return boost::math::cyl_bessel_j(fabs(order), x);
}

d_complex BesselJ::operator ()(const d_complex z)
{
    // Input values
    double zr = real(z);
    double zi = imag(z);
    double nu = fabs(order);
    int kode = 1;
    int N = 1;

    // Output values
    double cyr,cyi;
    int nz,ierr;

    // External function call
    zbesj_(&zr,&zi,&nu,&kode,&N,&cyr,&cyi,&nz,&ierr);
    d_complex answer(cyr,cyi);

    if (order < 0.0)
    {
        double cwrkr,cwrki,cyr2,cyi2;
        zbesy_(&zr,&zi,&nu,&kode,&N,&cyr2,&cyi2,&nz,&cwrkr,&cwrki,&ierr);

        d_complex answer2(cyr2,cyi2);

        return answer*cos(arma::datum::pi*nu)
                -answer2*sin(arma::datum::pi*nu);
    }

    return answer;
}
